% ptu_file_extract.m
%
% Written by James Thornham
%
% Script to call the read_ptu_v1 function from PicoQuant and extract
% individual channels, plot them over time, and create histograms.
% read_ptu_v1 only works for T3 Measure mode data, hence this program only

% works for T3 data.

% Section 1: Import data

clear
clc

% Import file and open using read_ptu_v1. Saves the data as "output".

[name,filepathlocation]=uigetfile('*.ptu');
totalfilepath=append(filepathlocation,name);
[output, txtout]=Read_PTU_V1(totalfilepath);

% IMPORTANT: You must change the sync cycle time. If the led is turned on
% every 750 us, the sync time will be 750e6 ps.

sync_cycle_time=20e9;

%% Section 2: Extract the individual channels (assumes only 0 and 1 are being used)

% This section creates a structure called "raw". The raw data is extracted
% from the output structure and values saved inside the raw structure can 
% be accessed using "raw.____". The following are saved inside raw:
% raw.tadj - The time from the start of the experiment at which each photon is detected
% raw.tadj_c1 - Adjusted times of photon arrival for channel 1 (raw.tadj_c0 is same for channel 0)
% raw.dtime_c1 - The time from the beginning of the last sync period, i.e. unadjusted time (raw.dtime_c0 is same for channel 0)
% raw.sync_c1 - sync number at which each photon occurs (raw.sync_c0 is same for channel 0)
% raw.binresolution - The time between data points in ps. The hydraharp creates bins with width of raw.binresolution and can detect 1 photon per bin
% raw.number_bins - The number of bins in each sync cycle (i.e. the maximum number of photons that can be detected per sync cycle)
% raw.Counts - Number of photons detected during the whole experiment for each channel


raw.tadj=(output.ph_dtime+(output.ph_sync-1)*(sync_cycle_time/output.Headers.MeasDesc_BinningFactor))*output.Headers.MeasDesc_BinningFactor;

[raw.Counts,]=groupcounts(output.ph_channel);

raw.tadj_c1=raw.tadj(output.ph_channel==1);
raw.tadj_c0=raw.tadj(output.ph_channel==0);

raw.dtime_c1=output.ph_dtime(output.ph_channel==1);
raw.dtime_c0=output.ph_dtime(output.ph_channel==0);
raw.sync_c1=output.ph_sync(output.ph_channel==1);
raw.sync_c0=output.ph_sync(output.ph_channel==0);
raw.binresolution=output.MeasDesc_Resolution*10^12;
raw.number_bins=sync_cycle_time/output.Headers.MeasDesc_BinningFactor;

% If you want to plot a circle at each time a photon is detected,
% un-comment the following line:

%plot_incidences(raw)

%% Section 3: Cut and bin the data

% A new structure called 'cut' is created where you can choose how much of
% each sync cycle to cut out. The values saved inside cut are the same as
% for raw except for there is no cut.tadj. There are other values added to
% cut that are not in raw:
% cut.sync_cycle_time - the time in ps per cycle. This is the same as sync_cycle_time which is set in Section 1.
% cut.cut_time - The time to remove from the start of each sync cycle in ms
% cut.binwidth - How many data points to use when plotting a histogram of all sync cycles
% cut.sum_c1 - The total counts over the summed sync cycles for channel 1 (cut.cum_c0 is the same for channel 0)
% cut.t_sum - The times associated with the cut.sum values

cut.cut_time=12.15e3; % Time to cut from the front of the data in usec
cut.binwidth=20;    % Binwidth to use for histograms. This is in number of 
                    % data points, not time. 
                    % To get time; bintime=binwidth*resolution 
                    % where resolution is set as the bin resolution during 
                    % the experiment.

cut.dtime_c1=raw.dtime_c1(raw.dtime_c1>=cut.cut_time/(raw.binresolution*10^(-6)));
cut.dtime_c0=raw.dtime_c0(raw.dtime_c0>=cut.cut_time/(raw.binresolution*10^(-6)));
cut.sync_c1=raw.sync_c1(raw.dtime_c1>=cut.cut_time/(raw.binresolution*10^(-6)));
cut.sync_c0=raw.sync_c0(raw.dtime_c0>=cut.cut_time/(raw.binresolution*10^(-6)));

cut.Counts=[size(cut.dtime_c0,1),size(cut.dtime_c1,1)];

cut.binresolution=raw.binresolution;
cut.number_bins=raw.number_bins;
cut.sync_cycle_time=sync_cycle_time;


% cut_dtime=output.ph_dtime(output.ph_dtime>(cut_time/0.065536));
% cut_channel=output.ph_channel(output.ph_dtime>(cut_time/0.065536));
% cut_sync=output.ph_sync(output.ph_dtime>(cut_time/0.065536));
% 
% [countspersync,syncnum]=groupcounts(output.ph_sync);
% [cut_countspersync,cut_syncnum]=groupcounts(cut_sync);

[times_0,times_1,Counts_c0,Counts_c1]=single_bin_hist(cut);

binned_hist(cut,cut.binwidth);

extract_lifetime(cut,cut.binwidth);

% To change the number of cycles to sum over, adjust the second value
% bellow. Remember, this is sync cycles, not time.

cut=total_counts(cut,10500);

%% Sum over a given number of sync cycles

function cut=total_counts(cut,cycles_to_sum)

    % Bins the number of photons measured in cycles_to_sum number of sync
    % cycles.

    [cut.sum_c1,cut.t_sum]=histcounts(cut.sync_c1,binwidth=cycles_to_sum);
    [cut.sum_c0,]=histcounts(cut.sync_c0,binwidth=cycles_to_sum);
    
    % Adjust cut.t_sum to be in ms
    cut.t_sum=cut.t_sum*(cut.binresolution*cut.number_bins/1e9);
    
    % Plots data
    figure
    
    plot(cut.t_sum(2:end-1),cut.sum_c1(1:end-1),'b')
    %hold on
    %plot(cut.t_sum(2:end-1),cut.sum_c0(1:end-1),'k')
    
    title(['Sum of counts over ',int2str(cycles_to_sum),' sync cycles'])
    xlabel('Time (ms)')
    ylabel('Counts')
    
    box off
    set(gca,TickDir='out')
    
    legend('Channel 1','Channel 0',EdgeColor='None')

end

%% Plot the Incidences over time for both channels

function plot_incidences(raw)
    
    % plots a + every time a photon arrives for both channels.
    figure
    
    plot(raw.tadj_c0,zeros(size(raw.tadj_c0,1),1),'k+',MarkerFaceColor='k',MarkerSize=1)
    hold on
    plot(raw.tadj_c1,ones(size(raw.tadj_c1,1),1),'b+',MarkerFaceColor='b',MarkerSize=1)
    
    title('Incidence Counts')
    xlabel('Time (ps)')
    ylabel('Channel')
    
    box off
    set(gca,TickDir='out')
    xlim([0 1000000000000])
    ylim([-0.1 1.1])
    
end

%% Plot histograms for each channel with bin size of 1 data point

function [times_0,times_1,Counts_c0,Counts_c1]=single_bin_hist(cut)
    
    % Find how many photons arrive at each time point during a sync cycle.
    [Counts_c1,times_1]=groupcounts(cut.dtime_c1);
    [Counts_c0,times_0]=groupcounts(cut.dtime_c0);
    
    % Convert times from data point to time in ps
    times_1=times_1*cut.binresolution;
    times_0=times_0*cut.binresolution;

    % Plots the histogram for both channels
    figure
    
    plot(times_1,Counts_c1,'b')
    hold on
    plot(times_0,Counts_c0,'k')
    
    title('Histogram (1 point bin)')
    xlabel('Time (ps)')
    ylabel('Counts')
    
    box off
    set(gca,TickDir='out')
    xlim([0 cut.sync_cycle_time])
    
    legend('Channel 1','Channel 0',EdgeColor='None')

end

%% Plot histogram with set bin size

function binned_hist(cut,binwidth)

    % Calculate the time width of each histogram bin
    bintime=binwidth*cut.binresolution;
    
    nbins=round(cut.number_bins/binwidth);
    
    % Chopping up the data. The edges array is 1 larger than the number of bins
    % since it lists 0 and the first value. The top end of each bin is used
    % as the x value to plot data.
    
    [C1Hist,edges_c1]=histcounts(cut.dtime_c1,nbins);
    [C0Hist,edges_c0]=histcounts(cut.dtime_c0,nbins);
    edges_c1=edges_c1*cut.binresolution;
    edges_c0=edges_c0*cut.binresolution;    

    figure
    
    plot(edges_c1(2:end),C1Hist,'b')
    hold on
    plot(edges_c0(2:end),C0Hist,'k')
    
    title(['Histogram: ',num2str(nbins),' point bins (',num2str(bintime/1000),' ns)'])
    xlabel('Time (ps)')
    ylabel('Counts')
    
    box off
    set(gca,TickDir='out')
    xlim([0 cut.sync_cycle_time])
    
    legend('Channel 1','Channel 0',EdgeColor='None')

end

%% Plot total photons per sync cycle against time

intensity_c0=accumarray(cut.sync_c0,1);
time=[sync_cycle_time/10^9:sync_cycle_time/10^9:size(intensity_c0,1)*sync_cycle_time/10^9]';
time=time/1000;
sum_intensity_c0=sum(intensity_c0) % sum of counts in measurement window channel 0

plot(time,intensity_c0);

intensity_c1=accumarray(cut.sync_c1,1);
sum_intensity_c1=sum(intensity_c1); % sum of counts in measurement window channel 1

%% Plot histograms for Individual Cycles
  % Plotting a histogram for oil cycle #1427
  % Select the cycle number from raw.sync_0
  oil_time=raw.dtime_c0(749486:749739); % Input the row numbers that correspond to the first and last event of that cycle number.
  oil_time=oil_time*raw.binresolution/1e9; % Convert x-axis from ps to ms
  oil_hist=histogram(oil_time,linspace(0,5768*raw.binresolution/1e9, 4769)); % Plot histogram of 4769 bins
  oil_edges=oil_hist.BinEdges; % Convert BinEdges in the histogram directory to a usable vector.
  oil_midpoints=(oil_edges(1:end-1)+oil_edges(2:end))/2; % Calculate bin midpoints
  oil_int=oil_hist.Values; % Retrieve y-values for bin height (photon counts)

  % Plotting a histogram for climb cycle #1489

  climb_time=raw.dtime_c0(765207:765735); % Input the row numbers that correspond to the first and last event of that cycle number.
  climb_time=climb_time*raw.binresolution/1e9; % Convert x-axis from ps to ms
  climb_hist=histogram(climb_time,linspace(0,5768*raw.binresolution/1e9, 4769)); % Plot histogram
  climb_edges=climb_hist.BinEdges; % Convert BinEdges in the histogram directory to a usable vector.
  climb_midpoints=(climb_edges(1:end-1)+climb_edges(2:end))/2; % Calculate bin midpoints
  climb_int=climb_hist.Values; % Retrieve y-values for bin height (photon counts)

  % Plotting a histogram for peak cycle #1524
  peak_time=raw.dtime_c0(765207:765735); % Input the row numbers that correspond to the first and last event of that cycle number.
  peak_time=peak_time*raw.binresolution/1e9; % Convert x-axis from ps to ms
  peak_hist=histogram(peak_time,linspace(0,5768*raw.binresolution/1e9, 4769)); % Plot histogram
  peak_edges=peak_hist.BinEdges; % Convert BinEdges in the histogram directory to a usable vector.
  peak_midpoints=(peak_edges(1:end-1)+peak_edges(2:end))/2; % Calculate bin midpoints
  peak_int=peak_hist.Values; % Retrieve y-values for bin height (photon counts)